# JCR REPLICATION for "How Civilian Loyalties Shape Rebel-led Victimization of Rebel Constituencies" #
# ILAYDA B. ONDER #
# MAIN ANALYSIS #
# MAY 27, 2024 #
# This .R file generates Table 4 and Figure 3.

load("Analysis Data.RData") # load analysis data

# TABLE 4 ----
#install.packages("rdd")
#install.packages("rddtools")
#install.packages("rdrobust")

library(rdd) # package for RDD analysis
library(rddtools) # package for RDD analysis
library(rdrobust) # package for RDD analysis

covs1 = "govrep_vil_log" # control variables for simpler models
covs2 = "govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + 
pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + 
y2014 + y2015 + y2016 + y2017 + y2018" # control variables for extended models

publicworker <- rdd_data(y = publicworker, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0) # create RDD data
collaborator <- rdd_data(y = collaborator, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0) # create RDD data
politician <- rdd_data(y = politician, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0) # create RDD data

# Run models
f1 <- rdd_reg_lm(rdd_object = publicworker, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covar.opt = list(strategy = "include")) # Model 1 of Table 4
se1 <- sqrt(diag(vcovHC(f1, type = "HC1"))) # store robust standard errors

f2 <- rdd_reg_lm(rdd_object = publicworker, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs1, covar.opt = list(strategy = "include")) # Model 2 of Table 4
se2 <- sqrt(diag(vcovHC(f2, type = "HC1"))) # store robust standard errors

f3 <- rdd_reg_lm(rdd_object = publicworker, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs2, covar.opt = list(strategy = "include")) # Model 3 of Table 4
se3 <- sqrt(diag(vcovHC(f3, type = "HC1"))) # store robust standard errors

f4 <- rdd_reg_lm(rdd_object = collaborator, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covar.opt = list(strategy = "include")) # Model 4 of Table 4
se4 <- sqrt(diag(vcovHC(f4, type = "HC1"))) # store robust standard errors

f5 <- rdd_reg_lm(rdd_object = collaborator, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covariates = covs1, covar.opt = list(strategy = "include")) # Model 5 of Table 4
se5 <- sqrt(diag(vcovHC(f5, type = "HC1"))) # store robust standard errors

f6 <- rdd_reg_lm(rdd_object = collaborator, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covariates = covs2, covar.opt = list(strategy = "include")) # Model 6 of Table 4
se6 <- sqrt(diag(vcovHC(f6, type = "HC1"))) # store robust standard errors

f7 <- rdd_reg_lm(rdd_object = politician, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covar.opt = list(strategy = "include")) # Model 7 of Table 4
se7 <- sqrt(diag(vcovHC(f7, type = "HC1"))) # store robust standard errors

f8 <- rdd_reg_lm(rdd_object = politician, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs1, covar.opt = list(strategy = "include")) # Model 8 of Table 4
se8 <- sqrt(diag(vcovHC(f8, type = "HC1"))) # store robust standard errors

f9 <- rdd_reg_lm(rdd_object = politician, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs2, covar.opt = list(strategy = "include")) # Model 9 of Table 4
se9 <- sqrt(diag(vcovHC(f9, type = "HC1"))) # store robust standard errors

# Export Table 4
stargazer::stargazer(f1, f2, f3, f4, f5, f6, f7, f8, f9, 
                     type = "text",
                     se = list(se1, se2, se3, se4, se5, se6, se7, se8, se9), 
                     keep = c("D", "govrep_vil_log", "govrep_non_log"),
                     covariate.labels = c("Treatment (Incumbent Party Victory)", "Government Violence", "Government Repression"),
                     column.labels = c(rep("Public Workers", 3), rep("Traitors", 3), rep("Pro-Gov Politicians", 3)),
                     dep.var.labels.include = F,
                     omit.stat = c("f", "ser"),
                     title = "Table 4: Support for the Government and the PKK’s Targeting of Specific Constituency Members, 2014-2019")




# FIGURE 3----
h1 = rddtools::rdd_bw_ik(kernel = "Normal", publicworker) # store I-K optimal bandwidth calculation
h2 = rddtools::rdd_bw_ik(kernel = "Normal", collaborator) # store I-K optimal bandwidth calculation
h3 = rddtools::rdd_bw_ik(kernel = "Normal", politician) # store I-K optimal bandwidth calculation

#install.packages("ggpubr")
#install.packages("ggthemr)")

library(ggpubr) # load package for graphs
library(ggthemr) # load package for graph themes
ggthemr('fresh') # set ggplot2 theme

covs2_fig = df_analysis[c('govrep_vil_log', 'govrep_non_log',
              'frequencyclash', 'intensityclash', 'murders', 'deaths', 'recruitswe', 'voters', 
              'pop', 'urbanpop', 'kurdvote', 'kurdpop', 'distance', 'altitude', 'rebellion', 'literacy', 'soceco',
              'y2014', 'y2015', 'y2016', 'y2017', 'y2018')] # control variables for extended models

# Figure 3 Panel A
fig3_a = rdplot(y = df_analysis$publicworker, x = df_analysis$margin_kurd, c = 0, p = 1, h = h1, covs = covs2_fig, title = "",
                x.lim = c(-h1, h1), y.lim = c(0, 0.15), x.label = "", y.label = "Outcome", covs_eval = "mean") 
fig3_a$rdplot$layers[[1]]$aes_params$size = 3 # graph specifications (e.g., size, color)
fig3_a$rdplot$layers[[2]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_a$rdplot$layers[[3]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_a$rdplot$layers[[4]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_a$rdplot$layers[[2]]$aes_params$colour = "firebrick" # graph specifications (e.g., size, color)
fig3_a$rdplot$layers[[3]]$aes_params$colour = "firebrick" # graph specifications (e.g., size, color)
fig3_a$rdplot + theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
                      plot.title = element_text(hjust = 0.5, face = "bold"),
                      axis.title = element_text(size = 18),
                      axis.text = element_text(size = 18),
                      legend.title = element_blank(),
                      legend.position = "bottom",
                      legend.key.height =  unit(0.5, 'cm'),
                      legend.key.width =  unit(4, 'cm'),
                      legend.text = element_text(size = 18)) # Generate Panel A of Figure 3

# Figure 3 Panel B
fig3_b = rdplot(y = df_analysis$collaborator, x = df_analysis$margin_kurd, c = 0, p = 1, h = h2,  covs = covs2_fig, title = "",
                x.lim = c(-h2, h2), y.lim = c(0, 0.15), x.label = "", y.label = "Outcome", covs_eval = "mean") 
fig3_b$rdplot$layers[[1]]$aes_params$size = 3 # graph specifications (e.g., size, color)
fig3_b$rdplot$layers[[2]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_b$rdplot$layers[[3]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_b$rdplot$layers[[4]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_b$rdplot$layers[[2]]$aes_params$colour = "firebrick" # graph specifications (e.g., size, color)
fig3_b$rdplot$layers[[3]]$aes_params$colour = "firebrick" # graph specifications (e.g., size, color)
fig3_b$rdplot + theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
                      plot.title = element_text(hjust = 0.5, face = "bold"),
                      axis.title = element_text(size = 18),
                      axis.text = element_text(size = 18),
                      legend.title = element_blank(),
                      legend.position = "bottom",
                      legend.key.height =  unit(0.5, 'cm'),
                      legend.key.width =  unit(4, 'cm'),
                      legend.text = element_text(size = 18)) # Generate Panel B of Figure 3

# Figure 3 Panel C
fig3_c = rdplot(y = df_analysis$politician, x = df_analysis$margin_kurd, c = 0, p = 1, h = h3, covs = covs2_fig, title = "",
            x.lim = c(-h3, h3), y.lim = c(0, 0.15), x.label = "", y.label = "Outcome")
fig3_c$rdplot$layers[[1]]$aes_params$size = 3 # graph specifications (e.g., size, color)
fig3_c$rdplot$layers[[2]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_c$rdplot$layers[[3]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_c$rdplot$layers[[4]]$aes_params$size = 2 # graph specifications (e.g., size, color)
fig3_c$rdplot$layers[[2]]$aes_params$colour = "firebrick" # graph specifications (e.g., size, color)
fig3_c$rdplot$layers[[3]]$aes_params$colour = "firebrick" # graph specifications (e.g., size, color)
fig3_c$rdplot + theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
                      plot.title = element_text(hjust = 0.5, face = "bold"),
                      axis.title = element_text(size = 18),
                      axis.text = element_text(size = 18),
                      legend.title = element_blank(),
                      legend.position = "bottom",
                      legend.key.height =  unit(0.5, 'cm'),
                      legend.key.width =  unit(4, 'cm'),
                      legend.text = element_text(size = 18))

# END ----




